#include "mt3d_cuda.h"

// Realize the virtual function of the solver class
void mt3d::AxProduct(cublasHandle_t cub_handle, cusparseHandle_t cus_handle, cusparseDnVecDescr_t x, 
    cusparseDnVecDescr_t prod_Ax, const int n_size, const int nz_size, cusparseOperation_t oper_t)
{
    // Calculate the product of A*x
    cusparseSpMV(cus_handle, oper_t, &_one, smat_A_, x, &_zero, prod_Ax, CUDA_C_64F, CUSPARSE_SPMV_ALG_DEFAULT, d_buf_);
    return;
}

void mt3d::MxProduct(cublasHandle_t cub_handle, cusparseHandle_t cus_handle, cusparseDnVecDescr_t x, 
    cusparseDnVecDescr_t prod_Ax, const int n_size, const int nz_size, cusparseOperation_t oper_t)
{
    // Preconditioning using triangular matrixes
    cusparseSpSV_solve(cus_handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &_one, smat_IC_, x, dvec_P_, 
        CUDA_C_64F, CUSPARSE_SPSV_ALG_DEFAULT, descr_L_);

    cusparseSpSV_solve(cus_handle, CUSPARSE_OPERATION_TRANSPOSE, &_one, smat_IC_, dvec_P_, prod_Ax, 
		CUDA_C_64F, CUSPARSE_SPSV_ALG_DEFAULT, descr_LT_);
    return;
}